package scu.maqiang.homogenization;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T101;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DoubleIntegerPair;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class SOTS3DT4101 {
	public SOTS3DT4101(FES3T41 fsT4, FES3T101 fsT10) {
		this.fsT4 = fsT4;
		this.fsT10 = fsT10;
	}
	
	/**
	 * ��������ȴ��������ͷ����ȴ���ϵ����һ�ס����׵������Լ����Ȼ��ȴ���ϵ����
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 * @param SOCS
	 */
	public void KxxCellSolution(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D, double[][] SOCS) {
		SRMatrix A = new SRMatrix(fsT4.GetNdof());
		fsT4.assembleStiff(kxxFunc, param, BVPType.COMMON, A);
		double[][] RHS = new double[3][fsT4.GetNdof()];
		fsT4.assembleGlobalVector(kxxFunc, param, this::elementkxxFirstOrder, BVPType.COMMON, RHS);
		fsT4.applyBC_MBN(A, Direct.All, 1, 2, 3, 4, 5, 6);
		fsT4.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
		fsT4.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
		NewIterSSolver solver = new NewIterSSolver(A);
		MVO.fill(FOCS, 0.0);
		solver.PCGSSOR(RHS[0], FOCS[0], 1.5, 1);
		solver.PCGSSOR(RHS[1], FOCS[1], 1.5, 1);
		solver.PCGSSOR(RHS[2], FOCS[2], 1.5, 1);
		
		computeHomogenizedKxx(kxxFunc, param, FOCS, D);
		
		RHS = new double[9][fsT4.GetNdof()];
		fsT4.assembleGlobalVector(kxxFunc, param, FOCS, D, this::elementKxxSecondOrder, BVPType.COMMON, RHS);
		for(int i = 0; i < 9; i++) {
			Arrays.fill(SOCS[i], 0.0);
			fsT4.applyBC_MBN(RHS[i], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
			solver.PCGSSOR(RHS[i], SOCS[i], 1.5, 1);
		}
	}
	
	/**
	 * �����ȴ���ϵ��һ�׵��������Ҷ���������
	 * @param xyz
	 * @param coef
	 * @param tp
	 * @param eleXYZ
	 */
	public void elementkxxFirstOrder(int[] ele, double[][] xyz, double[][] coef, BVPType tp, double[][] eleXYZ) {
		double[][] grad = new double[4][3];
		double volume = CG.gradient3DT4(xyz, grad);
		double temp = 0.25 * MVO.sum(coef[0]) * volume;

		for(int i = 0; i < 3; i++) {
			for(int j = 0; j < 4; j++) {
				eleXYZ[i][j] = -temp * grad[j][i];
			}
		}
//		eleXYZ[0][0] = -temp * grad[0][0];
//		eleXYZ[0][1] = -temp * grad[1][0];
//		eleXYZ[0][2] = -temp * grad[2][0];
//		eleXYZ[0][3] = -temp * grad[3][0];
//		eleXYZ[1][0] = -temp * grad[0][1];
//		eleXYZ[1][1] = -temp * grad[1][1];
//		eleXYZ[1][2] = -temp * grad[2][1];
//		eleXYZ[1][3] = -temp * grad[3][1];
//		eleXYZ[2][0] = -temp * grad[0][2];
//		eleXYZ[2][1] = -temp * grad[1][2];
//		eleXYZ[2][2] = -temp * grad[2][2];
//		eleXYZ[2][3] = -temp * grad[3][2];
	}
	
	/**
	 * �����ȴ���ϵ�����׵����������Ҷ���������
	 * @param xyz
	 * @param kxx
	 * @param FOCS
	 * @param homo
	 * @param tp
	 * @param eleXYZ
	 */
	public void elementKxxSecondOrder(int[] ele, double[][] xyz, double[][] kxx, double[][] FOCS, double[][] homo, BVPType tp, double[][] eleXYZ) {
		double[][] grad = new double[4][3];
		double volume = CG.gradient3DT4(xyz, grad);
		double[] dx = MVO.col(grad, 0);
		double[] dy = MVO.col(grad, 1);
		double[] dz = MVO.col(grad, 2);
		double N1Dx = MVO.dot_product(dx, FOCS[0]);
		double N1Dy = MVO.dot_product(dy, FOCS[0]);
		double N1Dz = MVO.dot_product(dz, FOCS[0]);
		double N2Dx = MVO.dot_product(dx, FOCS[1]);
		double N2Dy = MVO.dot_product(dy, FOCS[1]);
		double N2Dz = MVO.dot_product(dz, FOCS[1]);
		double N3Dx = MVO.dot_product(dx, FOCS[2]);
		double N3Dy = MVO.dot_product(dy, FOCS[2]);
		double N3Dz = MVO.dot_product(dz, FOCS[2]);
		double mKxx = MVO.average(kxx[0]);
		for(int i = 0; i < 4; i++) {
			eleXYZ[0][i] = (mKxx - homo[0][0] + mKxx * N1Dx - mKxx * MVO.sum(FOCS[0]) * dx[i]) * volume / 4.0; //11
			eleXYZ[1][i] = (     - homo[0][1] + mKxx * N1Dy - mKxx * MVO.sum(FOCS[0]) * dy[i]) * volume / 4.0; //12
			eleXYZ[2][i] = (     - homo[0][2] + mKxx * N1Dz - mKxx * MVO.sum(FOCS[0]) * dz[i]) * volume / 4.0; //13
			eleXYZ[3][i] = (     - homo[1][0] + mKxx * N2Dx - mKxx * MVO.sum(FOCS[1]) * dx[i]) * volume / 4.0; //21
			eleXYZ[4][i] = (mKxx - homo[1][1] + mKxx * N2Dy - mKxx * MVO.sum(FOCS[1]) * dy[i]) * volume / 4.0; //22
			eleXYZ[5][i] = (     - homo[1][2] + mKxx * N2Dz - mKxx * MVO.sum(FOCS[1]) * dz[i]) * volume / 4.0; //23
			eleXYZ[6][i] = (     - homo[2][0] + mKxx * N3Dx - mKxx * MVO.sum(FOCS[2]) * dx[i]) * volume / 4.0; //31
			eleXYZ[7][i] = (     - homo[2][1] + mKxx * N3Dy - mKxx * MVO.sum(FOCS[2]) * dy[i]) * volume / 4.0; //32
			eleXYZ[8][i] = (mKxx - homo[2][2] + mKxx * N3Dz - mKxx * MVO.sum(FOCS[2]) * dz[i]) * volume / 4.0; //33
		}
	}
	
	/**
	 * ��һ�׵������������ȴ������Ȼ�ϵ������
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 */
	public void computeHomogenizedKxx(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D) {
		double[][] grad = new double[4][3];
		double[][] xyz = new double[4][3];
		Mesh3T4 mesh = (Mesh3T4)fsT4.getMesh();
		double[] NN1 = new double[4];
		double[] NN2 = new double[4];
		double[] NN3 = new double[4];
		double N1Dx, N1Dy, N1Dz, N2Dx, N2Dy, N2Dz, N3Dx, N3Dy, N3Dz;
		double wholeVolume = 0.0;
		double mKxx = 0.0;
		double[][] coef = new double[kxxFunc.length][4];
		MVO.fill(D, 0.0);
		for(int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, xyz);
			MVO.getSlice(FOCS[0], ele, NN1);
			MVO.getSlice(FOCS[1], ele, NN2);
			MVO.getSlice(FOCS[2], ele, NN3);
			double volume = CG.gradient3DT4(xyz, grad);
			double[] dx = MVO.col(grad, 0);
			double[] dy = MVO.col(grad, 1);
			double[] dz = MVO.col(grad, 2);
			N1Dx = MVO.dot_product(NN1, dx);
			N1Dy = MVO.dot_product(NN1, dy);
			N1Dz = MVO.dot_product(NN1, dz);
			N2Dx = MVO.dot_product(NN2, dx);
			N2Dy = MVO.dot_product(NN2, dy);
			N2Dz = MVO.dot_product(NN2, dz);
			N3Dx = MVO.dot_product(NN3, dx);
			N3Dy = MVO.dot_product(NN3, dy);
			N3Dz = MVO.dot_product(NN3, dz);
			wholeVolume += volume;
			fsT4.getCoefFromFunc(kxxFunc, xyz, mesh.getElementLabel(i), param, coef);
			mKxx = MVO.average(coef[0]);
			D[0][0] += mKxx * (1 + N1Dx) * volume;
			D[0][1] += mKxx * N2Dx * volume;
			D[0][2] += mKxx * N3Dx * volume;
			D[1][0] += mKxx * N1Dy * volume;
			D[1][1] += mKxx * (1 + N2Dy) * volume;
			D[1][2] += mKxx * N3Dy * volume;
			D[2][0] += mKxx * N1Dz * volume;
			D[2][1] += mKxx * N2Dz * volume;
			D[2][2] += mKxx * (1 + N3Dz) * volume;
		}
		MVO.mul(D, 1.0 / wholeVolume);
	}
	
	public void assembleHeatSOTSSolution(FES3T41 fsFine, double[] T_0, double[] T_1, double[] T_2,
                                         FES3T101 fsHomo, double[] T_Homo, FES3T41 fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		double[] xy = new double[3];
		double[] cellXY = new double[3];
		
		double[] eleVal10 = new double[10];
		double[] eleVal4 = new double[4];

		double[][] grad = new double[3][fsHomo.GetNdof()];
		double[][] Se = new double[9][fsHomo.GetNdof()];
		//fsHomo.computeGradient(T_Homo, grad);
		
		fsHomo.computeFirstAndSecondDerivative(T_Homo, grad, Se);

		double corrector1, corrector2;
		
		DoubleIntegerPair homoData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdx10 = fsHomo.mesh.getElement(0).clone();
		
		DoubleIntegerPair cellData = new DoubleIntegerPair(0.0, 0);
		int[] eleIdx4 = fsCell.mesh.getElement(0).clone();
		double val = 0.0;
		for (int i = 0, n = fsFine.GetNdof(); i < n; i++) {
			fsFine.mesh.getCoordPerNode(i, xy);
			cellXY[0] = xy[0] * N - (int) (xy[0] * N);
			cellXY[1] = xy[1] * N - (int) (xy[1] * N);
			cellXY[2] = xy[2] * N - (int) (xy[2] * N);
			fsHomo.getValue(T_Homo, xy, homoData, eleIdx10, eleVal10);
			T_0[i] = homoData.first;
			fsCell.getValue(FOCS[0], cellXY, cellData, eleIdx4, eleVal4);
			val = cellData.first;
			
//			T_0[i] = fsHomo.getValue(T_Homo, xy, eleIdx6, eleVal6);
//			val = fsCell.getValue(FOCS[0], cellXY, eleIdx3, eleVal3);
			// corrector1 = 0.0;
			// corrector2 = 0.0;
			// corrector1 += val * MVO.dot_product(MVO.getSlice(grad[0], eleIdx6), eleVal6);
			// corrector1 += MVO.dot_product(MVO.getSlice(FOCS[1], eleIdx3), eleVal3)
			// * MVO.dot_product(MVO.getSlice(grad[1], eleIdx6), eleVal6);
			// T_1[i] = T_0[i] + ep * corrector1;
			// corrector2 += MVO.dot_product(MVO.getSlice(SOCS[0], eleIdx3), eleVal3) *
			// MVO.dot_product(MVO.getSlice(Se[0], eleIdx6), eleVal6)
			// + MVO.dot_product(MVO.getSlice(SOCS[1], eleIdx3), eleVal3) *
			// MVO.dot_product(MVO.getSlice(Se[1], eleIdx6), eleVal6)
			// + MVO.dot_product(MVO.getSlice(SOCS[2], eleIdx3), eleVal3) *
			// MVO.dot_product(MVO.getSlice(Se[2], eleIdx6), eleVal6)
			// + MVO.dot_product(MVO.getSlice(SOCS[3], eleIdx3), eleVal3) *
			// MVO.dot_product(MVO.getSlice(Se[3], eleIdx6), eleVal6);

			corrector1 = 0.0;
			corrector2 = 0.0;
			corrector1 += val * MVO.dot_product(grad[0], eleIdx10, eleVal10);
			corrector1 += MVO.dot_product(FOCS[1], eleIdx4, eleVal4) * MVO.dot_product(grad[1], eleIdx10, eleVal10);
			corrector1 += MVO.dot_product(FOCS[2], eleIdx4, eleVal4) * MVO.dot_product(grad[2], eleIdx10, eleVal10);
			T_1[i] = T_0[i] + ep * corrector1;
			corrector2 += MVO.dot_product(SOCS[0], eleIdx4, eleVal4) * MVO.dot_product(Se[0], eleIdx10, eleVal10);
			corrector2 += MVO.dot_product(SOCS[1], eleIdx4, eleVal4) * MVO.dot_product(Se[2], eleIdx10, eleVal10);
			corrector2 += MVO.dot_product(SOCS[2], eleIdx4, eleVal4) * MVO.dot_product(Se[2], eleIdx10, eleVal10);
			corrector2 += MVO.dot_product(SOCS[3], eleIdx4, eleVal4) * MVO.dot_product(Se[1], eleIdx10, eleVal10);
			corrector2 += MVO.dot_product(SOCS[4], eleIdx4, eleVal4) * MVO.dot_product(Se[0], eleIdx10, eleVal10);
			corrector2 += MVO.dot_product(SOCS[5], eleIdx4, eleVal4) * MVO.dot_product(Se[2], eleIdx10, eleVal10);
			corrector2 += MVO.dot_product(SOCS[6], eleIdx4, eleVal4) * MVO.dot_product(Se[2], eleIdx10, eleVal10);
			corrector2 += MVO.dot_product(SOCS[7], eleIdx4, eleVal4) * MVO.dot_product(Se[1], eleIdx10, eleVal10);
			corrector2 += MVO.dot_product(SOCS[8], eleIdx4, eleVal4) * MVO.dot_product(Se[0], eleIdx10, eleVal10);
			// Arrays.fill(cellValue, 0.0);
			// Arrays.fill(homoDxxValue, 0.0);
			// for(int j = 0; j < 3; j++) {
			// cellValue[0] += SOCS[0][eleIdx3[j]] * eleVal3[j];
			// cellValue[1] += SOCS[1][eleIdx3[j]] * eleVal3[j];
			// cellValue[2] += SOCS[2][eleIdx3[j]] * eleVal3[j];
			// cellValue[3] += SOCS[3][eleIdx3[j]] * eleVal3[j];
			// }
			// for(int j = 0; j < 6; j++) {
			// homoDxxValue[0] += Se[0][eleIdx6[j]] * eleVal6[j];
			// homoDxxValue[1] += Se[1][eleIdx6[j]] * eleVal6[j];
			// homoDxxValue[2] += Se[2][eleIdx6[j]] * eleVal6[j];
			// homoDxxValue[3] += Se[3][eleIdx6[j]] * eleVal6[j];
			// }
			T_2[i] = T_1[i] + ep * ep * corrector2;
		}
	}

	FES3T41 fsT4;
	FES3T101 fsT10;
}
